/* CS580 Homework 3 */

#include	"stdafx.h"
#include	"stdio.h"
#include	"math.h"
#include	"Gz.h"
#include	"rend.h"


int GzRotXMat(float degree, GzMatrix mat)
{
// Create rotate matrix : rotate along x axis
// Pass back the matrix using mat value
	mat[1][1]=cos(degree*PI/180);
	mat[1][2]=-sin(degree*PI/180);
	mat[2][1]=sin(degree*PI/180);
	mat[2][2]=cos(degree*PI/180);

	return GZ_SUCCESS;
}


int GzRotYMat(float degree, GzMatrix mat)
{
// Create rotate matrix : rotate along y axis
// Pass back the matrix using mat value
	mat[0][0]=cos(degree*PI/180);
	mat[2][0]=-sin(degree*PI/180);
	mat[0][2]=sin(degree*PI/180);
	mat[2][2]=cos(degree*PI/180);

	return GZ_SUCCESS;
}


int GzRotZMat(float degree, GzMatrix mat)
{
// Create rotate matrix : rotate along z axis
// Pass back the matrix using mat value
	mat[0][0]=cos(degree*PI/180);
	mat[0][1]=-sin(degree*PI/180);
	mat[1][0]=sin(degree*PI/180);
	mat[1][1]=cos(degree*PI/180);

	return GZ_SUCCESS;
}


int GzTrxMat(GzCoord translate, GzMatrix mat)
{
// Create translation matrix
// Pass back the matrix using mat value
	mat[0][3]=translate[0];
	mat[1][3]=translate[1];
	mat[2][3]=translate[2];

	return GZ_SUCCESS;
}


int GzScaleMat(GzCoord scale, GzMatrix mat)
{
// Create scaling matrix
// Pass back the matrix using mat value
	mat[0][0]=scale[0];
	mat[1][1]=scale[1];
	mat[2][2]=scale[2];

	return GZ_SUCCESS;
}


//----------------------------------------------------------
// Begin main functions

int GzNewRender(GzRender **render, GzRenderClass renderClass, GzDisplay	*display)//
{
/*  
- malloc a renderer struct /
- keep closed until all inits are done /
- setup Xsp and anything only done once /
- span interpolator needs pointer to display /
- check for legal class GZ_Z_BUFFER_RENDER /
- init default camera /
*/ 
	int i,j,k;
	GzCoord x,y,z;
	float cl,norm;

	(*render) = (GzRender*) operator new (sizeof(GzRender));
	(*render)->renderClass = renderClass;
	(*render)->display = display;
	(*render)->ambientlight = NULL;
	(*render)->interp_mode = 0;
	(*render)->matlevel = 0;
	(*render)->numlights = 1;
	(*render)->open = display->open;
	(*render)->spec = DEFAULT_SPEC;
	(*render)->tex_fun = NULL;

	(*render)->Ka[RED] = 0.1;			/*throw DEFAULT_AMBIENT, DEFAULT_DIFFUSE, DEFAULT_SPECULAR to Ka, Kd, Ks*/
	(*render)->Ka[GREEN] = 0.1;
	(*render)->Ka[BLUE] = 0.1;
	(*render)->Kd[RED] = 0.7;
	(*render)->Kd[GREEN] = 0.6;
	(*render)->Kd[BLUE] = 0.5;
	(*render)->Ks[RED] = 0.2;
	(*render)->Ks[GREEN] = 0.3;
	(*render)->Ks[BLUE] = 0.4;

	/*Initialize camera by defualt*/
	(*render)->camera.FOV = DEFAULT_FOV;
	(*render)->camera.lookat[X] = 0.0;				/*Init lookat position*/
	(*render)->camera.lookat[Y] = 0.0;
	(*render)->camera.lookat[Z] = 0.0;
	(*render)->camera.position[X] = DEFAULT_IM_X;	/*Init world position*/
	(*render)->camera.position[Y] = DEFAULT_IM_Y;
	(*render)->camera.position[Z] = DEFAULT_IM_Z;
	(*render)->camera.worldup[X] = 0.0;				/*Init world up vector*/
	(*render)->camera.worldup[Y] = 1.0;
	(*render)->camera.worldup[Z] = 0.0;
	for(i=0;i<4;i++)								/*Init Xpi matrix*/
	{
		for(j=0;j<4;j++)
		{
			if(i==j)
			{
				(*render)->camera.Xpi[i][j] = 1.0;
			}
			else
			{
				(*render)->camera.Xpi[i][j] = 0.0;
			}
		}
	}
	(*render)->camera.Xpi[3][2] = tan(((*render)->camera.FOV/2)*PI/180);

	for(i=0;i<4;i++)								/*Init Xiw matrix*/
	{
		for(j=0;j<4;j++)
		{
			if(i==j)
			{
				(*render)->camera.Xiw[i][j] = 1.0;
			}
			else
			{
				(*render)->camera.Xiw[i][j] = 0.0;
			}
		}
	}

	/*According to prof, which is SideVec = ViewDir x WorldUp, and UpVector = SideVec x ViewDir*/
	cl = sqrt(((*render)->camera.lookat[X]-(*render)->camera.position[X])*((*render)->camera.lookat[X]-(*render)->camera.position[X])
					+((*render)->camera.lookat[Y]-(*render)->camera.position[Y])*((*render)->camera.lookat[Y]-(*render)->camera.position[Y])
					+((*render)->camera.lookat[Z]-(*render)->camera.position[Z])*((*render)->camera.lookat[Z]-(*render)->camera.position[Z]));
	z[X] = ((*render)->camera.lookat[X]-(*render)->camera.position[X])/cl;
	z[Y] = ((*render)->camera.lookat[Y]-(*render)->camera.position[Y])/cl;
	z[Z] = ((*render)->camera.lookat[Z]-(*render)->camera.position[Z])/cl;				/*Init Z*/
	norm = sqrt((z[Z]*(*render)->camera.worldup[Y]-z[Y]*(*render)->camera.worldup[Z])*(z[Z]*(*render)->camera.worldup[Y]-z[Y]*(*render)->camera.worldup[Z])
				+(z[X]*(*render)->camera.worldup[Z]-z[Z]*(*render)->camera.worldup[X])*(z[X]*(*render)->camera.worldup[Z]-z[Z]*(*render)->camera.worldup[X])
				+(z[Y]*(*render)->camera.worldup[X]-z[X]*(*render)->camera.worldup[Y])*(z[Y]*(*render)->camera.worldup[X]-z[X]*(*render)->camera.worldup[Y]));
	x[X] = -(z[Y]*(*render)->camera.worldup[Z]-z[Z]*(*render)->camera.worldup[Y])/norm;
	x[Y] = -(z[Z]*(*render)->camera.worldup[X]-z[X]*(*render)->camera.worldup[Z])/norm;
	x[Z] = -(z[X]*(*render)->camera.worldup[Y]-z[Y]*(*render)->camera.worldup[X])/norm;	/*Init X*/
	norm = sqrt((x[Z]*z[Y]-x[Y]*z[Z])*(x[Z]*z[Y]-x[Y]*z[Z])+(x[X]*z[Z]-x[Z]*z[X])*(x[X]*z[Z]-x[Z]*z[X])+(x[Y]*z[X]-x[X]*z[Y])*(x[Y]*z[X]-x[X]*z[Y]));
	y[X] = -(x[Y]*z[Z]-x[Z]*z[Y])/norm;
	y[Y] = -(x[Z]*z[X]-x[X]*z[Z])/norm;
	y[Z] = -(x[X]*z[Y]-x[Y]*z[X])/norm;													/*Init Y*/
	/*******************************************************************************************/

	/*Compute Xiw*/
	(*render)->camera.Xiw[0][0] = x[X];
	(*render)->camera.Xiw[0][1] = x[Y];
	(*render)->camera.Xiw[0][2] = x[Z];
	(*render)->camera.Xiw[1][0] = y[X];
	(*render)->camera.Xiw[1][1] = y[Y];
	(*render)->camera.Xiw[1][2] = y[Z];
	(*render)->camera.Xiw[2][0] = z[X];
	(*render)->camera.Xiw[2][1] = z[Y];
	(*render)->camera.Xiw[2][2] = z[Z];
	(*render)->camera.Xiw[0][3] = -(x[X]*(*render)->camera.position[X]+x[Y]*(*render)->camera.position[Y]+x[Z]*(*render)->camera.position[Z]);
	(*render)->camera.Xiw[1][3] = -(y[X]*(*render)->camera.position[X]+y[Y]*(*render)->camera.position[Y]+y[Z]*(*render)->camera.position[Z]);
	(*render)->camera.Xiw[2][3] = -(z[X]*(*render)->camera.position[X]+z[Y]*(*render)->camera.position[Y]+z[Z]*(*render)->camera.position[Z]);
	/******************************/

	/*Initialize color*/
	for(i=0;i<3;i++)
	{
		(*render)->flatcolor[i] = 0;
	}

	/*Initialize the lights*/
	for(i=0;i<MAX_LIGHTS;i++)
	{
		(*render)->lights[i] = NULL;
	}

	/*Initialize matrix stack*/
	for(i=0;i<MATLEVELS;i++)
	{
		for(j=0;j<4;j++)
		{
			for(k=0;k<4;k++)
			{
				if(k==j)
				{
					(*render)->Ximage[i][j][k] = 1.0;
					(*render)->Xnorm[i][j][k] = 1.0;
				}
				else
				{
					(*render)->Ximage[i][j][k] = 0.0;
					(*render)->Xnorm[i][j][k] = 0.0;
				}
			}
		}
	}

	/*Initialize Xsp matrix*/
	for(i=0;i<4;i++)
	{
		for(j=0;j<4;j++)
		{
			(*render)->Xsp[i][j] = 0.0;
		}
	}
	(*render)->Xsp[0][0] = display->xres/2;
	(*render)->Xsp[0][3] = display->xres/2;
	(*render)->Xsp[1][1] = -display->yres/2;
	(*render)->Xsp[1][3] = display->yres/2;
	(*render)->Xsp[2][2] = INT_MAX*tan(((*render)->camera.FOV/2)*PI/180);
	(*render)->Xsp[3][3] = 1.0;
	/***********************/

	if(!GZ_Z_BUFFER_RENDER)
	{
		return GZ_FAILURE;
	}
	return GZ_SUCCESS;

}


int GzFreeRender(GzRender *render)
{
/* 
-free all renderer resources
*/
//	GzFreeDisplay(render->display);
	delete(render);
	return GZ_SUCCESS;
}


int GzBeginRender(GzRender *render)
{
/*  
- set up for start of each frame - clear frame buffer /
- compute Xiw and projection xform Xpi from camera definition /
- init Ximage - put Xsp at base of stack, push on Xpi and Xiw /
- now stack contains Xsw and app can push model Xforms if it want to. /
*/
	int i,j,k;
	GzCoord x,y,z;
	float cl,norm;

	GzInitDisplay(render->display);					/*Set up for start of each frame*/

	/*Compute Xiw and Xpi*/
	for(i=0;i<4;i++)								/*Init Xpi matrix*/
	{
		for(j=0;j<4;j++)
		{
			if(i==j)
			{
				render->camera.Xpi[i][j] = 1.0;
			}
			else
			{
				render->camera.Xpi[i][j] = 0.0;
			}
		}
	}
	render->camera.Xpi[3][2] = tan((render->camera.FOV/2)*PI/180);

	for(i=0;i<4;i++)								/*Init Xiw matrix*/
	{
		for(j=0;j<4;j++)
		{
			if(i==j)
			{
				render->camera.Xiw[i][j] = 1.0;
			}
			else
			{
				render->camera.Xiw[i][j] = 0.0;
			}
		}
	}
	/*According to prof, which is SideVec = ViewDir x WorldUp, and UpVector = SideVec x ViewDir*/
	cl = sqrt((render->camera.lookat[X]-render->camera.position[X])*(render->camera.lookat[X]-render->camera.position[X])
					+(render->camera.lookat[Y]-render->camera.position[Y])*(render->camera.lookat[Y]-render->camera.position[Y])
					+(render->camera.lookat[Z]-render->camera.position[Z])*(render->camera.lookat[Z]-render->camera.position[Z]));
	z[X] = (render->camera.lookat[X]-render->camera.position[X])/cl;
	z[Y] = (render->camera.lookat[Y]-render->camera.position[Y])/cl;
	z[Z] = (render->camera.lookat[Z]-render->camera.position[Z])/cl;				/*Init Z*/
	norm = sqrt((z[Z]*render->camera.worldup[Y]-z[Y]*render->camera.worldup[Z])*(z[Z]*render->camera.worldup[Y]-z[Y]*render->camera.worldup[Z])
				+(z[X]*render->camera.worldup[Z]-z[Z]*render->camera.worldup[X])*(z[X]*render->camera.worldup[Z]-z[Z]*render->camera.worldup[X])
				+(z[Y]*render->camera.worldup[X]-z[X]*render->camera.worldup[Y])*(z[Y]*render->camera.worldup[X]-z[X]*render->camera.worldup[Y]));
	x[X] = -(z[Y]*render->camera.worldup[Z]-z[Z]*render->camera.worldup[Y])/norm;
	x[Y] = -(z[Z]*render->camera.worldup[X]-z[X]*render->camera.worldup[Z])/norm;
	x[Z] = -(z[X]*render->camera.worldup[Y]-z[Y]*render->camera.worldup[X])/norm;	/*Init X*/
	norm = sqrt((x[Z]*z[Y]-x[Y]*z[Z])*(x[Z]*z[Y]-x[Y]*z[Z])+(x[X]*z[Z]-x[Z]*z[X])*(x[X]*z[Z]-x[Z]*z[X])+(x[Y]*z[X]-x[X]*z[Y])*(x[Y]*z[X]-x[X]*z[Y]));
	y[X] = -(x[Y]*z[Z]-x[Z]*z[Y])/norm;
	y[Y] = -(x[Z]*z[X]-x[X]*z[Z])/norm;
	y[Z] = -(x[X]*z[Y]-x[Y]*z[X])/norm;												/*Init Y*/
	/*******************************************************************************************/

	/*Compute Xiw*/
	render->camera.Xiw[0][0] = x[X];
	render->camera.Xiw[0][1] = x[Y];
	render->camera.Xiw[0][2] = x[Z];
	render->camera.Xiw[1][0] = y[X];
	render->camera.Xiw[1][1] = y[Y];
	render->camera.Xiw[1][2] = y[Z];
	render->camera.Xiw[2][0] = z[X];
	render->camera.Xiw[2][1] = z[Y];
	render->camera.Xiw[2][2] = z[Z];
	render->camera.Xiw[0][3] = -(x[X]*render->camera.position[X]+x[Y]*render->camera.position[Y]+x[Z]*render->camera.position[Z]);
	render->camera.Xiw[1][3] = -(y[X]*render->camera.position[X]+y[Y]*render->camera.position[Y]+y[Z]*render->camera.position[Z]);
	render->camera.Xiw[2][3] = -(z[X]*render->camera.position[X]+z[Y]*render->camera.position[Y]+z[Z]*render->camera.position[Z]);
	/*********************/

	/*Update Xsp matrix for change of FOV*/
	render->Xsp[2][2] = INT_MAX*tan((render->camera.FOV/2)*PI/180);
	/*************************************/

	/*Initialize Ximage matrix stack*/
	for(i=0;i<MATLEVELS;i++)
	{
		for(j=0;j<4;j++)
		{
			for(k=0;k<4;k++)
			{
				if(k==j)
				{
					render->Ximage[i][j][k] = 1.0;
				}
				else
				{
					render->Ximage[i][j][k] = 0.0;
				}
			}
		}
	}
	/********************************/

	/*Initialize Ximage matrix stack*/
	GzPushMatrix(render,render->Xsp);			/*Push Xsp*/
	GzPushMatrix(render,render->camera.Xpi);	/*Push Xpi*/
	GzPushMatrix(render,render->camera.Xiw);	/*Push Xiw*/
	/********************************/

	return GZ_SUCCESS;
}

int GzPutCamera(GzRender *render, GzCamera *camera)
{
/*
- overwrite renderer camera structure with new camera definition
*/
	render->camera.FOV = camera->FOV;					/*Overwrite FOV*/
	render->camera.lookat[X] = camera->lookat[X];		/*Overwrite lookat position*/
	render->camera.lookat[Y] = camera->lookat[Y];
	render->camera.lookat[Z] = camera->lookat[Z];
	render->camera.position[X] = camera->position[X];	/*Overwrite world position*/
	render->camera.position[Y] = camera->position[Y];
	render->camera.position[Z] = camera->position[Z];
	render->camera.worldup[X] = camera->worldup[X];		/*Overwrite world up vector*/
	render->camera.worldup[Y] = camera->worldup[Y];
	render->camera.worldup[Z] = camera->worldup[Z];

	return GZ_SUCCESS;	
}

int GzPushMatrix(GzRender *render, GzMatrix	matrix)
{
/*
- push a matrix onto the Ximage stack
- check for stack overflow
*/
	int i,j;
	if(render->matlevel==0)
	{
		for(i=0;i<4;i++)
		{
			for(j=0;j<4;j++)
			{
				render->Ximage[render->matlevel][i][j]=matrix[i][j];
			}
		}
	}
	else
	{
		MatrixMultiply(render->Ximage[render->matlevel-1],matrix,render->Ximage[render->matlevel]);
	}
	render->matlevel++;
	if(render->matlevel>=MATLEVELS)
	{
		return GZ_FAILURE;
	}

	return GZ_SUCCESS;
}

int GzPopMatrix(GzRender *render)
{
/*
- pop a matrix off the Ximage stack
- check for stack underflow
*/
	/*Matlevel minus one first, then pop out the matrix*/
	render->matlevel--;
	if(render->matlevel<0)
	{
		return GZ_FAILURE;
	}

	return GZ_SUCCESS;
}


int GzPutAttribute(GzRender	*render, int numAttributes, GzToken	*nameList, GzPointer *valueList) /* void** valuelist */
{
/*
- set renderer attribute states (e.g.: GZ_RGB_COLOR default color)
- later set shaders, interpolaters, texture maps, and lights
*/
	for(int i=0;i<numAttributes;i++)
	{
		if(nameList[i]==GZ_RGB_COLOR)//set color
		{
			float* color=(float*)valueList[i];
			render->flatcolor[RED]=color[RED];
			render->flatcolor[GREEN]=color[GREEN];
			render->flatcolor[BLUE]=color[BLUE];
		}
	}
	
	return GZ_SUCCESS;
}

int GzPutTriangle(GzRender	*render, int numParts, GzToken *nameList, GzPointer *valueList)//
/* numParts : how many names and values */
{
/*  
- pass in a triangle description with tokens and values corresponding to 
      GZ_POSITION:3 vert positions in model space /
- Xform positions of verts  /
- Clip - just discard any triangle with verts behind view plane /
       - test for triangles with all three verts off-screen /
- invoke triangle rasterizer  /
*/
	float delta;
	double a,b,c,d;
	double z;
	int i;
	float tmp[4][1]={1.0,1.0,1.0,1.0};
	float result[4][1]={1.0,1.0,1.0,1.0};

	//HW2 w/ Xform and Clip
	for(i=0;i<numParts;i++)
	{
		if(nameList[i]==GZ_NULL_TOKEN)
		{
		}
		else if(nameList[i]==GZ_POSITION)//set triangle
		{
			GzCoord vertex[3];
			float* temp=(float*)valueList[i];
			GzIntensity oldr,oldg,oldb,olda;
			GzDepth oldZ;
			//initial edge DDAs and span DDAs
			EdgeDDA* e1=(EdgeDDA*) operator new (sizeof(EdgeDDA));
			EdgeDDA* e2=(EdgeDDA*) operator new (sizeof(EdgeDDA));
			EdgeDDA* e3=(EdgeDDA*) operator new (sizeof(EdgeDDA));
			SpanDDA* s=(SpanDDA*) operator new (sizeof(SpanDDA));
			//get the location of each point
			vertex[0][X]=temp[0];
			vertex[0][Y]=temp[1];
			vertex[0][Z]=temp[2];
			vertex[1][X]=temp[3];
			vertex[1][Y]=temp[4];
			vertex[1][Z]=temp[5];
			vertex[2][X]=temp[6];
			vertex[2][Y]=temp[7];
			vertex[2][Z]=temp[8];

			/*Need to apply matrix here*/
			for(i=0;i<3;i++)
			{
				tmp[0][0]=vertex[i][X];
				tmp[1][0]=vertex[i][Y];
				tmp[2][0]=vertex[i][Z];
				MatrixVectorMultiply(render->Ximage[render->matlevel-1],tmp,result);
				vertex[i][X]=result[0][0]/result[3][0];
				vertex[i][Y]=result[1][0]/result[3][0];
				vertex[i][Z]=result[2][0]/result[3][0];
			}
			/***************************/

			/*Discard any triangle with verts behind view plane*/
			if((vertex[0][Z]<0)||(vertex[1][Z]<0)||(vertex[2][Z]<0))
			{
				continue;
			}
			/***************************************************/

			//get the color
			short red = ctoi(render->flatcolor[RED]);
			short green = ctoi(render->flatcolor[GREEN]);
			short blue = ctoi(render->flatcolor[BLUE]);

			/*Rasterization start*/
			YSort(e1,e2,e3,vertex);//sort verts by Y
			LRTest(e1,e2,e3);//setup edge DDAs and sort them by L or R
			//calculate (a,b,c) of ax+by+cz+d=0
			/*calculate a,b,c for the triangle*/
			a=(int)((((double)vertex[1][Y]-(double)vertex[0][Y])*((double)vertex[2][Z]-(double)vertex[0][Z]))-(((double)vertex[1][Z]-(double)vertex[0][Z])*((double)vertex[2][Y]-(double)vertex[0][Y])));
			b=(int)((((double)vertex[1][Z]-(double)vertex[0][Z])*((double)vertex[2][X]-(double)vertex[0][X]))-(((double)vertex[1][X]-(double)vertex[0][X])*((double)vertex[2][Z]-(double)vertex[0][Z])));
			c=(int)((((double)vertex[1][X]-(double)vertex[0][X])*((double)vertex[2][Y]-(double)vertex[0][Y]))-(((double)vertex[1][Y]-(double)vertex[0][Y])*((double)vertex[2][X]-(double)vertex[0][X])));
			/**/

			if((e1->lr!=2)&&(e2->lr!=3))
			{
				AdvanceYCeiling(e1,e3);//Advance and DDA current positions to top y-scan line (ceiling)

				while(e1->vCurrent[Y]<=e1->vEnd[Y])
				{//start from the lower edge
					InitSpanDDA(s,e1,e3);//Initial span DDA
					d=((a*s->xCurrent+b*e1->vCurrent[Y]+c*s->zCurrent)+(a*s->xEnd+b*e1->vCurrent[Y]+c*s->zEnd))/2;//calculate d

					delta=(int)s->xStart+1-s->xStart;
					if(delta!=1)
					{
						s->xCurrent=s->xStart+delta;
						s->zCurrent=s->zStart+s->slopez*delta;
					}
					else
					{
						s->xCurrent=s->xStart;
						s->zCurrent=s->zStart;
					}//setup span edges and set positions

					while(s->xCurrent<=s->xEnd)
					{
						z=(d-a*(double)s->xCurrent-b*(double)e1->vCurrent[Y])/c;//interpolate z
						s->zCurrent=z;
						GzGetDisplay(render->display,s->xCurrent,e1->vCurrent[Y],&oldr,&oldg,&oldb,&olda,&oldZ);//get the z-buffer
						if(z<=oldZ)//check the z with the z-buffer
						{
							GzPutDisplay(render->display,(int)s->xCurrent,(int)e1->vCurrent[Y],red,green,blue,1,z);//test interpolate-Z and write color value into FB
						}
						s->xCurrent+=1;//next pixel
						s->zCurrent+=s->slopez;
					}
					//next line
					e1->vCurrent[Y]+=1;
					e1->vCurrent[X]=e1->vCurrent[X]+e1->slopex;
					e1->vCurrent[Z]=e1->vCurrent[Z]+e1->slopez;
					e3->vCurrent[Y]+=1;
					e3->vCurrent[X]=e3->vCurrent[X]+e3->slopex;
					e3->vCurrent[Z]=e3->vCurrent[Z]+e3->slopez;
				}

				delta=(int)e2->vStart[Y]+1-e2->vStart[Y];
				if(delta!=1)
				{
					e2->vCurrent[X]=e2->vCurrent[X]+e2->slopex*delta;
					e2->vCurrent[Y]=e2->vCurrent[Y]+delta;
					e2->vCurrent[Z]=e2->vCurrent[Z]+e2->slopez*delta;
				}//go to higher egde

				while(e3->vCurrent[Y]<=e3->vEnd[Y])
				{//start the scan line between the higher edge and the othe edge
					InitSpanDDA(s,e2,e3);//Initial span DDA
					d=((a*s->xCurrent+b*e3->vCurrent[Y]+c*s->zCurrent)+(a*s->xEnd+b*e3->vCurrent[Y]+c*s->zEnd))/2;//calculate d
					delta=(int)s->xStart+1-s->xStart;
					if(delta!=1)
					{
						s->xCurrent=s->xStart+delta;
						s->zCurrent=s->zStart+s->slopez*delta;
					}
					else
					{
						s->xCurrent=s->xStart;
						s->zCurrent=s->zStart;
					}//setup span edges and set positions

					while(s->xCurrent<=s->xEnd)
					{
						z=(d-a*(double)s->xCurrent-b*(double)e2->vCurrent[Y])/c;//interpolate z
						s->zCurrent=z;
						GzGetDisplay(render->display,s->xCurrent,e2->vCurrent[Y],&oldr,&oldg,&oldb,&olda,&oldZ);//get the z-buffer
						if(z<=oldZ)//check the z with the z-buffer
						{
							GzPutDisplay(render->display,(int)s->xCurrent,(int)e2->vCurrent[Y],red,green,blue,1,z);//test interpolate-Z and write color value into FB
						}
						s->xCurrent+=1;//next pixel
						s->zCurrent+=s->slopez;
					}
					//next line
					e2->vCurrent[Y]+=1;
					e2->vCurrent[X]=e2->vCurrent[X]+e2->slopex;
					e2->vCurrent[Z]=e2->vCurrent[Z]+e2->slopez;
					e3->vCurrent[Y]+=1;
					e3->vCurrent[X]=e3->vCurrent[X]+e3->slopex;
					e3->vCurrent[Z]=e3->vCurrent[Z]+e3->slopez;
				}
			}
			else if(e1->lr==2)
			{
				AdvanceYCeiling(e2,e3);//Advance and DDA current positions to top y-scan line (ceiling)
				while(e3->vCurrent[Y]<=e3->vEnd[Y])
				{//since the top edge is horizontal, only the span edge between other two edges needs to be calculated
					InitSpanDDA(s,e3,e2);//Initial span DDA
					d=((a*s->xCurrent+b*e3->vCurrent[Y]+c*s->zCurrent)+(a*s->xEnd+b*e3->vCurrent[Y]+c*s->zEnd))/2;//calculate d
					delta=(int)s->xStart+1-s->xStart;
					if(delta!=1)
					{
						s->xCurrent=s->xStart+delta;
						s->zCurrent=s->zStart+s->slopez*delta;
					}
					else
					{
						s->xCurrent=s->xStart;
						s->zCurrent=s->zStart;
					}//setup span edges and set positions

					while(s->xCurrent<=s->xEnd)
					{
						z=(d-a*(double)s->xCurrent-b*(double)e2->vCurrent[Y])/c;//interpolate z
						s->zCurrent=z;
						GzGetDisplay(render->display,s->xCurrent,e2->vCurrent[Y],&oldr,&oldg,&oldb,&olda,&oldZ);//get the z-buffer
						if(z<=oldZ)//check the z with the z-buffer
						{
							GzPutDisplay(render->display,(int)s->xCurrent,(int)e2->vCurrent[Y],red,green,blue,1,z);//test interpolate-Z and write color value into FB
						}
						s->xCurrent+=1;//next pixel
						s->zCurrent+=s->slopez;
					}
					//next line
					e2->vCurrent[Y]+=1;
					e2->vCurrent[X]=e2->vCurrent[X]+e2->slopex;
					e2->vCurrent[Z]=e2->vCurrent[Z]+e2->slopez;
					e3->vCurrent[Y]+=1;
					e3->vCurrent[X]=e3->vCurrent[X]+e3->slopex;
					e3->vCurrent[Z]=e3->vCurrent[Z]+e3->slopez;
				}
			}
			else if(e2->lr==3)
			{
				AdvanceYCeiling(e1,e3);//Advance and DDA current positions to top y-scan line (ceiling)
				while(e1->vCurrent[Y]<=e1->vEnd[Y])
				{//since the bottom edge is horizontal, only the span edge between other two edges needs to be calculated
					InitSpanDDA(s,e1,e3);//Initial span DDA
					d=((a*s->xCurrent+b*e1->vCurrent[Y]+c*s->zCurrent)+(a*s->xEnd+b*e1->vCurrent[Y]+c*s->zEnd))/2;//calculate d
					delta=(int)s->xStart+1-s->xStart;
					if(delta!=1)
					{
						s->xCurrent=s->xStart+delta;
						s->zCurrent=s->zStart+s->slopez*delta;
					}
					else
					{
						s->xCurrent=s->xStart;
						s->zCurrent=s->zStart;
					}//setup span edges and set positions

					while(s->xCurrent<=s->xEnd)
					{
						z=(d-a*s->xCurrent-b*e1->vCurrent[Y])/c;//interpolate z
						s->zCurrent=z;
						GzGetDisplay(render->display,s->xCurrent,e1->vCurrent[Y],&oldr,&oldg,&oldb,&olda,&oldZ);//get the z-buffer
						if(z<=oldZ)//check the z with the z-buffer
						{
							GzPutDisplay(render->display,(int)s->xCurrent,(int)e1->vCurrent[Y],red,green,blue,1,z);//test interpolate-Z and write color value into FB
						}
						s->xCurrent+=1;//next pixel
						s->zCurrent+=s->slopez;
					}
					//next line
					e1->vCurrent[Y]+=1;
					e1->vCurrent[X]=e1->vCurrent[X]+e1->slopex;
					e1->vCurrent[Z]=e1->vCurrent[Z]+e1->slopez;
					e3->vCurrent[Y]+=1;
					e3->vCurrent[X]=e3->vCurrent[X]+e3->slopex;
					e3->vCurrent[Z]=e3->vCurrent[Z]+e3->slopez;
				}
			}
			/*Rasterization end*/
			delete(e1);
			delete(e2);
			delete(e3);
			delete(s);
		}
	}

	return GZ_SUCCESS;
}

/* NOT part of API - just for general assistance */

short	ctoi(float color)		/* convert float color to GzIntensity short */
{
  return(short)((int)(color * ((1 << 12) - 1)));
}

void YSort(EdgeDDA *e1,EdgeDDA *e2,EdgeDDA *e3,GzCoord *vertex)
{//Sort verts by Y from low to high
	GzCoord temp;
	temp[X]=vertex[0][X];
	temp[Y]=vertex[0][Y];
	temp[Z]=vertex[0][Z];
	for(int i=0;i<3;i++)
	{
		for(int j=i;j<3;j++)
		{
			if(vertex[j][Y]<vertex[i][Y])
			{
				temp[X]=vertex[i][X];
				temp[Y]=vertex[i][Y];
				temp[Z]=vertex[i][Z];
				vertex[i][X]=vertex[j][X];
				vertex[i][Y]=vertex[j][Y];
				vertex[i][Z]=vertex[j][Z];
				vertex[j][X]=temp[X];
				vertex[j][Y]=temp[Y];
				vertex[j][Z]=temp[Z];
			}
			else if(vertex[j][Y]==vertex[i][Y])
			{
				if(vertex[j][X]<=vertex[i][X])
				{
					temp[X]=vertex[i][X];
					temp[Y]=vertex[i][Y];
					temp[Z]=vertex[i][Z];
					vertex[i][X]=vertex[j][X];
					vertex[i][Y]=vertex[j][Y];
					vertex[i][Z]=vertex[j][Z];
					vertex[j][X]=temp[X];
					vertex[j][Y]=temp[Y];
					vertex[j][Z]=temp[Z];
				}
			}
		}
	}
	//Initial edge DDAs
	if(vertex[0][Y]==vertex[1][Y])
	{//if the top edge is horizontal
		InitEdgeDDAT(e1,vertex[X],vertex[Y]);
		InitEdgeDDA(e2,vertex[Y],vertex[Z]);
		InitEdgeDDA(e3,vertex[X],vertex[Z]);

	}
	else if(vertex[1][Y]==vertex[2][Y])
	{//if the bottom edge is horizontal
		InitEdgeDDAB(e2,vertex[Y],vertex[Z]);
		InitEdgeDDA(e1,vertex[X],vertex[Y]);
		InitEdgeDDA(e3,vertex[X],vertex[Z]);
	}
	else
	{
		InitEdgeDDA(e1,vertex[X],vertex[Y]);
		InitEdgeDDA(e2,vertex[Y],vertex[Z]);
		InitEdgeDDA(e3,vertex[X],vertex[Z]);
	}
}

void InitEdgeDDA(EdgeDDA *e,float *v1,float *v2)
{//Initial edge DDA
	e->vStart[X]=v1[X];
	e->vStart[Y]=v1[Y];
	e->vStart[Z]=v1[Z];
	e->vEnd[X]=v2[X];
	e->vEnd[Y]=v2[Y];
	e->vEnd[Z]=v2[Z];
	e->vCurrent[X]=v1[X];
	e->vCurrent[Y]=v1[Y];
	e->vCurrent[Z]=v1[Z];
	e->slopex=(v2[X]-v1[X])/(v2[Y]-v1[Y]);
	e->slopez=(v2[Z]-v1[Z])/(v2[Y]-v1[Y]);
	e->lr=L;
}

void InitEdgeDDAT(EdgeDDA *e,float *v1,float *v2)
{//Initial horizontal top edge DDA 
	e->vStart[X]=v1[X];
	e->vStart[Y]=v1[Y];
	e->vStart[Z]=v1[Z];
	e->vEnd[X]=v2[X];
	e->vEnd[Y]=v2[Y];
	e->vEnd[Z]=v2[Z];
	e->vCurrent[X]=v1[X];
	e->vCurrent[Y]=v1[Y];
	e->vCurrent[Z]=v1[Z];
	e->slopex=0;
	e->slopez=0;
	e->lr=T;
}

void InitEdgeDDAB(EdgeDDA *e,float *v1,float *v2)
{//Initial horizontal bottom edge DDA
	e->vStart[X]=v1[X];
	e->vStart[Y]=v1[Y];
	e->vStart[Z]=v1[Z];
	e->vEnd[X]=v2[X];
	e->vEnd[Y]=v2[Y];
	e->vEnd[Z]=v2[Z];
	e->vCurrent[X]=v1[X];
	e->vCurrent[Y]=v1[Y];
	e->vCurrent[Z]=v1[Z];
	e->slopex=0;
	e->slopez=0;
	e->lr=B;
}

void LRTest(EdgeDDA *e1,EdgeDDA *e2,EdgeDDA *e3)
{//Sort edge DDA by L or R
	if(e1->lr==2)
	{//if top edge is horizontal
		e3->lr=L;
		e2->lr=R;
	}
	else if(e2->lr==3)
	{//if bottom edge is horizontal
		e1->lr=L;
		e3->lr=R;
	}
	else
	{
		if(e1->slopex<e3->slopex)
		{
			e3->lr=R;
		}
		else
		{
			e1->lr=R;
			e2->lr=R;
		}//sort L/R by slope
	}
}

void InitSpanDDA(SpanDDA *s,EdgeDDA *e1,EdgeDDA *e2)
{//Initial span DDA
	if(e1->lr==0)
	{//if edge1 is the left edge
		if(e1->vCurrent[X]<0)
		{//if the left end is out of the bound, start the scan line from x=0
			s->xStart=0;
			s->slopez=(e2->vCurrent[Z]-e1->vCurrent[Z])/(e2->vCurrent[X]-e1->vCurrent[X]);
			s->zStart=e1->vCurrent[Z]+s->slopez*(0-e1->vCurrent[X]);
			s->xEnd=e2->vCurrent[X];
			s->zEnd=e2->vCurrent[Z];
			s->xCurrent=s->xStart;
			s->zCurrent=s->zStart;
		}
		else
		{
			s->xStart=e1->vCurrent[X];
			s->zStart=e1->vCurrent[Z];
			s->xEnd=e2->vCurrent[X];
			s->zEnd=e2->vCurrent[Z];
			s->xCurrent=e1->vCurrent[X];
			s->zCurrent=e1->vCurrent[Z];
			s->slopez=(s->zEnd-s->zStart)/(s->xEnd-s->xStart);
		}
	}
	else if(e2->lr==0)
	{//if edge2 is left edge
		if(e2->vCurrent[X]<0)
		{//if the left end is out of the bound, start the scan line from x=0
			s->xStart=0;
			s->slopez=(e1->vCurrent[Z]-e2->vCurrent[Z])/(e1->vCurrent[X]-e2->vCurrent[X]);
			s->zStart=e2->vCurrent[Z]+s->slopez*(0-e2->vCurrent[X]);
			s->xEnd=e1->vCurrent[X];
			s->zEnd=e1->vCurrent[Z];
			s->xCurrent=s->xStart;
			s->zCurrent=s->zStart;
		}
		else
		{
			s->xStart=e2->vCurrent[X];
			s->zStart=e2->vCurrent[Z];
			s->xEnd=e1->vCurrent[X];
			s->zEnd=e1->vCurrent[Z];
			s->xCurrent=e2->vCurrent[X];
			s->zCurrent=e2->vCurrent[Z];
			s->slopez=(s->zEnd-s->zStart)/(s->xEnd-s->xStart);
		}
	}
}

void AdvanceYCeiling(EdgeDDA *e1,EdgeDDA *e2)
{//Advance and DDA current positions to top y-scan line (ceiling)
	float delta;
	delta=(int)e1->vStart[Y]+1-e1->vStart[Y];
	if(delta!=1)
	{
		e1->vCurrent[X]=e1->vCurrent[X]+e1->slopex*delta;
		e1->vCurrent[Y]=e1->vCurrent[Y]+delta;
		e1->vCurrent[Z]=e1->vCurrent[Z]+e1->slopez*delta;
	}
	delta=(int)e2->vStart[1]+1-e2->vStart[1];
	if(delta!=1)
	{
		e2->vCurrent[X]=e2->vCurrent[X]+e2->slopex*delta;
		e2->vCurrent[Y]=e2->vCurrent[Y]+delta;
		e2->vCurrent[Z]=e2->vCurrent[Z]+e2->slopez*delta;
	}
}

void MatrixMultiply (float (*a)[4], float (*b)[4], float (*c)[4])
{//Multiplication between two 4x4 matrixes, matrix a and b are the source and matrix c is the target
	int i,j,k;
	float tmp;

	for(i=0;i<4;i++)
	{
		for(j=0;j<4;j++)
		{
			tmp=0.0;
			for(k=0;k<4;k++)
			{
				tmp+=a[i][k]*b[k][j];
			}
			c[i][j]=tmp;
		}
	}
	return;
}

void MatrixVectorMultiply (float (*a)[4], float (*b)[1], float (*c)[1])
{//Multiplication between one 4x4 matrix and one 4x1 matrix, and matrix a is 4x4, matrix b is 4x1, and the target matrix c is 4x1
	int i,j,k;
	float tmp;

	for(i=0;i<4;i++)
	{
		for(j=0;j<1;j++)
		{
			tmp=0.0;
			for(k=0;k<4;k++)
			{
				tmp+=a[i][k]*b[k][j];
			}
			c[i][j]=tmp;
		}
	}
	return;
}

void CopyMatrix(float (*target)[4],float (*source)[4])
{//Copy one 4x4 matrix into another 4x4 matrix
	int i,j;

	for(i=0;i<4;i++)
	{
		for(j=0;j<4;j++)
		{
			target[i][j]=source[i][j];
		}
	}
}